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ABSTRACT 


Felix  Bloch's  1928  article  made  a  prediction  concerning  the  dynamical  behavior 
of  electrons  in  a  solid,  subject  to  a  uniform,  static  electric  field.  This  aspect  of  his  work, 
as  later  clarified  by  Zener,  showed  that  electrons  accelerated  by  an  electric  field  in  a 
periodic  potential,  under  the  right  conditions,  would  oscillate.  A  theoretical  debate  as  to 
the  existence  of  this  phenomenon  has  been  ongoing  since  Bloch's  proposal.  One  of  the 
most  controversial  consequences  of  this  prediction  is  that  an  electron  undergoing  Bloch 
oscillations  would  radiate.  The  controversy  on  the  theoretical  analysis  was  due  to  the 
great  difficulty  in  systematically  and  reliably  treating  interband  transitions  by  analytical 
methods  based  on  the  time-dependent  Schrodinger  equation  for  independent  electrons.  In 
this  thesis,  we  numerically  solve  the  time-dependent  Schrodinger  equation  to  show  that 
electrons  accelerated  by  in  an  electric  field  in  periodic  structures  do  undergo  Bloch 
oscillations  and  other  dynamic  behavior.  By  accurately  modeling  this  phenomenon,  we 
hope  to  gain  a  better  understanding  of  it  in  hopes  of  using  it  in  future  applications  as  a 
stable  source  of  Terahertz  (THz)  radiation. 
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I.  INTRODUCTION 


In  1928,  Felix  Bloch  [Ref.  1]  wrote  a  groundbreaking  work  discussing  the 
quantum  mechanical  nature  of  electrons  in  solids.  In  this  seminal  work,  Bloch 
demonstrated  what  is  widely  known  as  "Bloch's  theorem,"  which  establishes  the 
mathematical  form  of  the  electron  wave  function  of  a  crystalline  solid.  We  will  discuss 
Bloch's  theorem  later  in  this  thesis.  Bloch's  theorem  is  the  basis  for  our  understanding  of 
the  all-important  energy  bands  in  solids.  In  addition  to  establishing  Bloch's  theorem,  a 
less  well  known  result  from  Bloch's  1928  article  is  a  prediction  concerning  the  dynamical 
behavior  of  electrons  in  a  solid,  subject  to  a  uniform,  static  electric  field.  This  aspect  of 
his  work,  as  later  expanded  by  Zener  [Ref.  2],  showed  that  electrons  accelerated  by  an 
electric  field  in  a  periodic  potential,  under  the  right  conditions  would  oscillate.  The 
frequency  of  oscillation  is  given  by  the  formula  co  B  =  eFa/h  ,  where  e  is  the  charge  of  an 
electron,  F  is  the  electric  field  amplitude,  a  is  the  lattice  repeat  distance,  and  h  is 
Planck's  constant  h  divided  by  2:t .  This  phenomenon  has  come  to  be  known  as  Bloch 
oscillations.  A  theoretical  debate  as  to  the  existence  of  this  phenomenon  has  been 
ongoing  since  Bloch's  proposal.  One  of  the  most  controversial  consequences  of  this 
prediction  is  that  an  electron  undergoing  Bloch  oscillations  would  radiate.  The 
controversy  on  the  theoretical  analysis  was  due  to  the  great  difficulty  in  systematically 
and  reliably  treating  interband  transitions  by  analytical  methods  based  on  the  time- 
dependent  Schrodinger  equation  for  independent  electrons. 

With  the  advent  and  refinement  of  semiconductor  technology,  it  has  been  finally 
possible  to  test  the  predictions  of  Bloch's  proposal.  In  1970,  Esaki  and  Tsu  [Ref.  3] 
proposed  using  a  semiconductor  superlattice,  which  would  provide  the  needed  periodic 
structure  on  a  large  enough  scale  in  order  to  conduct  tests  to  generate  and  detect  Bloch 
oscillations.  According  to  Esaki  and  Tsu  [Ref.  3],  the  oscillating  electron  in  a 
semiconductor  superlattice  would  provide  a  source  of  Terahertz  (THz)  radiation.  It  was 
only  in  1992,  however,  that  Bloch  oscillations  and  the  associated  Terahertz  radiation  were 
conclusively  detected  in  semiconductor  superlattices  [Ref.  4].  These  recent  four  wave 


mixing  experiments  confirm  the  presence  of  photons  emitted  as  a  result  of  the  oscillating 
electron,  but  noted  the  oscillations  died  between  1  and  15  oscillations  [Ref.  4].  More 
recently,  Luban  and  co-workers[Ref.  5]  modeled  the  effect  of  an  electron  in  a 
semiconductor  superlattice  using  computer  techniques,  and  showed  Bloch  oscillations 
exist  at  Terahertz  frequencies,  but  did  not  show  the  decay  as  seen  in  real  world 
experiments.  The  advantage  of  a  purely  numerical  approach  to  solving  the  time- 
dependent  Schrodinger  equation  is  that  one  explicitly  works  with  the  full  Hamiltonian, 
and  not  some  approximation  as  in  analytical  solutions,  thus  including  all  interband 
transitions  in  the  simulation. 

This  paper  will  cover  the  theoretical  debates  and  arguments  surrounding  the  issue. 
Also,  it  will  develop  the  theory  behind  Bloch  oscillations.  Lastly,  it  will  cover  the 
numerical  techniques  involved  in  solving  the  Time  Dependent  Schrodinger  Equation,  and 
present  the  results  of  the  simulations. 


II.  THEORETICAL  WORK 

In  this  chapter,  we  summarize  and  discuss  the  works  of  major  contributors  on  the 
subject  of  Bloch  oscillations. 

A.  FELIX  BLOCH 

In  1928,  Felix  Bloch  wrote  a  ground  breaking  paper  [Ref.  1]  which  until  recently 
was  the  subject  of  great  debate.  In  it,  he  used  quantum  mechanics  to  model  the  effects  of 
electrons  in  periodic  (atomic)  structures.  In  this  landmark  paper,  he  proved  what  is  now 
known  as  Bloch's  theorem,  which  shows  that  the  mathematical  form  of  the  eigenvector  of 
the  time  independent  Schrodinger  equation  must  be: 

where 

KAx)  =  unAx  +  a)  i22) 

and  n  labels  the  energy  bands  and  is  known  as  the  band  index.  Also,  the  parameter  k  is 
the  allowed  wavevector  of  the  electron  and  is  restricted  to  the  Brillouin  zone  associated 
with  the  particular  lattice  structure  of  the  crystal.  The  associated  energy  eigenvalue 
En(k)  is  called  the  band  structure  function. 

Bloch  showed  [Ref.  1]  that  the  wave  function  of  any  particle  in  a  periodic 
potential  can  be  represented  as  a  sum  of  orthogonal  functions  that  have  the  periodicity  of 
the  lattice  structure.  Bloch  [Ref.  1],  starting  with  the  time  dependent  Schrodinger 
equation  with  an  external  electric  field  applied: 


,                2m                                    ilm  dxb(x,t) 
V2xp(x,t)-  —  (V(x)-eFx)V(x,t)  +  - ^-^  =  0  ,  (2.3) 


and  assuming  the  following  form  of  the  wave  function: 

-— E  i 

Mx,t)=2cnk(t)^ntke   h   '  ,  (2.4) 

njc 

derived  the  result  for  the  time  variation  of  the  amplitude  of  the  coefficients  of  the  Bloch 
functions  of  the  electron  wave  function: 

j;\cJ  =  -^f^\c,M-  (2-5) 

at  h     dk 

He  also  showed  for  a  collection  of  electrons  obeying  Fermi  statistics  [Ref.  1],  that  the 
time  variation  of  the  electron  distribution  function  f(k,t)  is: 

#_2raFtf  (2.6) 

dt  h     dk 

In  (2.6),  f(k,t)  is  the  probability  of  finding  an  electron  between  states  k  and  k  +  M  at 

time  t .    Bloch  claims  this  is  the  same  result  obtained  by  Lorentz  [Ref.  6]  in  an  earlier 
paper  published  in  1916. 

B.  CLARENCE  ZENER 

Zener  [Ref.  2],  in  1933,  substantially  clarified  Bloch's  work  [Ref.  1]  in  a  paper 
investigating  the  breakdown  of  solids  by  an  electric  field.    The  breakdown  of  a  dielectric 
under  the  influence  of  an  electric  field  takes  place  by  one  of  two  mechanisms  [Ref.  2]: 

1.  A  process  similar  to  the  electrical  breakdown  of  a  gas,  i.e.,  avalanche 
breakdown. 

2.  A  process  like  the  auto  ionization  of  free  atoms  in  an  electric  field. 

For  solids,  Zener  [Ref .  2]  said  the  second  process  is  more  important  in  a  solid  than  in  a 
gas  since  the  mean  free  path  of  an  electron  is  much  less  for  a  solid  than  a  gas.  An 


electron,  therefore,  tends  to  become  excited  and  move  about  in  the  lattice  instead  of  being 
freed  from  the  solid  completely.  He  used  the  Bloch  model  of  an  electron  to  calculate  the 
rate  at  which  electrons  move  from  lower  to  upper  energy  bands  in  a  solid  under  the 
influence  of  a  constant  electric  field. 

Zener  [Ref.  2]  started  with  the  same  initial  wavefunction  (2.1).  He  said  that  the 
condition  that  \i>k  be  finite  over  the  lattice  imposes  the  restriction  that  k  be  real  as  well. 

This  only  occurs  when  the  electron  energy  En  lies  within  the  allowed  energy  bands  that 
are  a  solution  of  the  time  independent  Schrodinger  equation,  with  the  first  energy  band 
occuring  within  the  values  -  n/a  s  k  <  n/a  [Ref.  2]  for  a  one  dimensional  model. 

Following  Bloch,  Zener  [Ref.  2]  expanded  his  wave  function  in  terms  of  wave  functions 
of  the  lowest  energy  band  with  no  electric  field  applied: 

ip(jt,0-  )ag(k,t)x^nk(x)dk.  (2.7) 

Following  Bloch's  derivation,  Zener  [Ref.  2]  reproduced  equation  (2.5)  in  a  more  concise 
form: 


j-M-^iM2  (2-8) 

at  h     dk 

where 


\g?-G[k-^ft),  (2.9) 

with  G  being  an  arbitrary  function.  Zener  [Ref.  2]  also  showed  that  the  velocity  in 
k  space  is  2xeF/h .  Using  the  limits  of  k ,  Zener  showed  the  period  of  time  for  real 
space  oscillation  is: 


h 

tb  =  (23t/a)/(2xeF/h)  = (2.10) 

eFa 


or  the  frequency  of  oscillation  is  given  by: 

vB=  —  .  (2.11) 

h 

Zener  [Ref.  2]  said  than  an  electron  confined  to  a  single  energy  band  will  move  opposite 
to  the  direction  of  the  field  until  reflected  by  the  lattice,  and  then  it  moves  in  the  opposite 
direction  until  it  is  stopped  by  the  field,  where  it  starts  the  same  motion  over  again. 
Zener  [Ref.  2]  then  calculated  the  probability  that  an  electron  would  undergo  transitions 
between  energy  bands.  Starting  with  a  wave  function  that  is  periodic  in  time,  Zener 
showed  the  wave  equation  for  an  electron  when  an  external  field  F  is  present  is  (2.3). 
Zener  [Ref.  2]  found  an  approximate  solution,  assuming  the  term  eFx  varies  slightly  over 
the  lattice  distance  a  : 


^,*(*)  =  ^M*)  (2-12) 


where  Kis  a  function  of  (E  -  eFx) .  Using  (2.12),  and  solving  foiK,  Zener  [Ref.  2] 
showed  the  rate  of  transition  (y  )  between  energy  bands  is  given  by: 


eFa 

y  = exp 


jt2  maE2^ 
~¥  \eF\  , 


(2.13) 


where  Eg  is  the  energy  gap  between  bands. 

C.  W.  V.  HOUSTON 

In  1939,  Houston  [Ref.  7]  showed  that  the  wave  vector  of  an  electron  accelerated 
by  an  electric  field  in  a  periodic  potential  grows  linearly  with  time  within  the  bounds  of 
the  first  Brillouin  zone[Ref.  7]. 


Houston  assumed  the  form  of  this  wave  function  as: 


-£fEk+\t<** 


Mx,t)  =  uk+^x)e'{k^e  (2.14) 


where  k  =  eF/h.  [Ref.  7],  which  shows  the  solution  of  the  time  dependent  Schrodinger 
equation  is  constructed  from  the  solution  of  the  time  independent  Schrodinger  equation, 
modulated  by  a  time  dependent  wave  vector.    Equation  (2.14)  is  a  solution  of  (2.3),  with 
the  following  remainder: 


jKWu^e^e  (2.15) 


The  remainder  term  is  zero  for  the  case  when  no  electric  field  is  present  (k  =  0),  and  is 

zero  when  Vwi+X,  =  0  for  the  free  electron  case  [Ref.  7].  In  the  case  where  X  or  Vm^  are 

not  zero,  but  their  product  is  very  small,  Houston  [Ref.  7]  claimed  (2.14)  is  a  good 
approximate  solution  to  (2.3).  In  the  case  where  the  remainder  is  not  small,  Houston 
[Ref.  7]  showed  that  the  electron  wave  function  can  be  written  as  an  infinite  sum  of 
orthogonal  Bloch  eigenfunctions: 


W*>0-2fl/(0"*+™-+Xr(*)*(  e  (2-16) 


where  /  are  reciprocal  lattice  vectors  associated  with  the  crystal  structure. 


Using  the  method  of  variation  of  constants,  Houston  solved  for  the  time  variation  of  the 
coefficients  a^t) .  Having  solved  this,  he  showed  [Ref.  7]  that  the  probability  for  the 

electron  to  transition  across  a  Brillouin  zone  boundary  is: 


P  =  4jize-a/(l  +  e-a)  (2.17) 

where  a  =  mV2a2/eFh2 .  Houston  claims  [Ref.  7]  this  result  is  similar  to  Zener's  [Ref. 
2]- 

D.  GREGORY  WANNIER 

Wannier's  contribution  to  the  field  of  solid  state  physics  is  immense.  One  of  his 
many  contributions  was  his  work  on  the  motion  of  electrons  in  solids  in  the  presence  of  a 
uniform  electric  field,  and  the  nature  of  the  energy  bands  formed  by  the  electric  field  in 
that  solid.    In  a  series  of  papers  from  1960  to  1968  [Refs.  8-11],  Wannier  decribed  the 
effects  of  electric  and  magnetic  .fields  on  electrons  in  periodic  potentials.  Wannier's  [Ref 
.8]  most  significant  contribution  was  to  show,  for  a  single  band  Hamiltonian,  that  the 
energy  levels  of  a  solid  under  the  effect  of  an  electric  field  become  discretized  such  that 
they  are  a  function  of  the  lattice  repeat  distance  a  : 

r./a 

E„=±  JWq(k)dk  +  nFa  (2.18) 

o 

where  [Ref.  9]  En  are  the  energy  levels,  the  integral  term  is  the  mean  energy  of  the  band, 

and  nFa  are  the  ladder  spacing  intervals.  This  series  of  equally  spaced  energy  levels  is 
termed  the  Stark  ladder,  or  sometimes  the  Wannier-Stark  ladder.  The  basic  idea  behind 
the  Stark  ladder  is  that  the  electric  field  destroys  the  degeneracy  of  the  atomic  energies  of 
the  solid  that  form  the  bands.  Thus,  in  Wannier's  picture,  electric  field  destroys  the  zero- 
field  energy  bands  and  replaces  them  with  a  ladder  of  equally  spaced  energy  levels.  The 
existence  of  Stark  ladders  in  solids  has  itself  been  as  controversial  a  topic  as  Bloch 


oscillations,  and  Stark  ladders  were  first  observed  in  semiconductor  superlattices  for  a 
single  band  beginning  in  the  1980's  [Ref.  12].    Wannier  also  showed  [Ref.  10]  that  if  the 
electron  wave  function  is  comprised  of  Bloch  functions,  then  the  energy  bands  are 
decoupled,  which  supports  his  idea  of  a  discrete  energy  bands. 

Perhaps  Wannier's  most  widely  known  contribution  to  solid  state  physics  are  the 
"Wannier  functions,"  which  are  in  essence  a  Fourier  transform  of  the  Bloch  states  [Ref. 
13].  Whereas  Bloch  states  are  extended  states,  i.e.,  the  electron  can  be  found  over  the 
entire  crystal,  the  Wannier  functions  are  localized  to  within  a  few  lattice  sites.  It  can  be 
shown  that  Wannier  functions  centered  about  different  lattice  sites  are  orthogonal  [Ref. 
13].  Also,  it  was  proved  that  Wannier  functions  at  the  same  lattice  site,  but  representing 
different  energy  bands  are  orthogonal  [Ref.  13].  The  Wannier  functions  are  thus  a 
convenient  basis  set  with  which  to  represent  the  electron  wave  function. 

E.  J.ZAK 

Zak  has  written  a  great  deal  on  the  subject  of  electron  motion  in  periodic 
potentials  under  the  influence  of  external  electric  and  magnetic  fields.  Probably  his  most 
outstanding  contribution  to  this  field  is  the  derivation  and  use  of  a  Bloch-type  set  of 
functions  which  are  expressible  as  localized  Wannier  functions  [Ref.  14].  The 
representation  of  these  functions  is  known  as  the  kq  representation  [Ref.  15],  and  was 

derived  using  the  fact  that  finite  vector  translations  in  real  and  reciprocal  space  in  a  lattice 
form  a  complete  set  of  commuting  observables.    The  difference  between  his  set  of 
functions  and  normal  Bloch-type  functions  used  to  solve  solid  state  problems  is  that  the 
former  is  a  function  of  two  continuous  variables  k  and  q ,  while  the  latter  are  function  of 

a  discrete  band  index  n  and  continuous  wave  vector  k  [Ref.  15].  He  was  an  outspoken 
critic,  however,  of  the  idea  of  a  Stark  ladder  in  a  solid  [Ref.  16].  His  main  criticism  was 
that  Wannier's  derivation  was  an  approximation  using  perturbation  on  a  one  band 
Hamiltonian  [Ref.  8],  and,  using  the  exact  time  independent  Schrodinger  equation,  Zak 
showed  the  integral  term  of  (2.19)  was,  in  fact,  an  arbitrary  constant,  thus  invalidating  the 


idea  of  a  ladder  of  energy  levels.  Wannier's  reply  [Ref.  17]  was  to  show  that  the  integral 
term  was  valid,  and  that  Stark  ladders  were  a  valid  model.  Zak's  reply  to  Wannier  [Ref. 
18]  was  that  Wannier  still  uses  an  approximate  method,  and  his  kq  approach  is  correct. 
Finally,  Zak  [Ref.  19]  claims,  with  a  more  robust  proof,  that  the  energy  spectrum  for  a 
Bloch  electron  in  an  electric  field  is  continuous.  The  author,  thus  far,  has  not  been  able  to 
determine  which  approach  is  correct.  Multi-band  Stark  levels  are  the  electrical  analog  to 
Landau  levels  in  a  magnetic  field,  and  ultimately  must  be  proved  or  disproved  by 
experimentalists. 

F.  A.  RABINOVITCH 

Rabinovitch's  first  comment  [Ref.  20]  on  the  effect  of  an  electric  field  on  an 
electron  in  a  solid  concerned  the  translational  symmetry  argument  in  finite  crystals. 
Rabinovitch  [Ref.  20]  showed  Born-von  Karman  periodic  boundary  conditions  for  finite 
crystals  in  an  electric  field  are  incompatible  with  the  Schrodinger  equation;  that  is,  they 
do  not  lead  to  the  expected  Stark  ladder  result  for  the  energy  levels.  He  [Ref.  20]  also 
showed  any  boundary  condition  compatible  with  the  Schrodinger  equation  breaks 
translational  symmetry,  and  prevents  the  use  of  Bloch  electron  functions.  His  other 
works  are  collaborations  with  Zak  [Refs.  21-22].  In  [Ref.  22],  using  Zak's  kq 
representation  for  the  Schrodinger  equation,  and  using  an  identical  single  band 
approximation  as  Houston,  Rabinovitch  and  Zak  were  able  to  derive  the  same  result  as 
Houston  (2.15)  for  the  electron  wave  function.  They  argued,  since  this  single  band 
approximation  ignores  other  bands  and  interband  coupling,  Houston's  solution  is 
unphysical  and  meaningless.  Also,  they  proposed  that  an  electron  composed  of  the  wave 
function  in  (2.15)  would  only  oscillate  with  time  faster  than  Inh/eFa ,  and  thus  Bloch 
oscillations  cannot  exist. 
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G.  J.  N.  CHURCHILL  AND  F.  E.  HOLMSTROM 

In  a  series  of  three  papers  [Refs.  23-25],  Churchill  and  Holmstrom  comment  on 
the  subject  of  Bloch  oscillations  and  electron  motion  in  periodic  potentials  subject  to  a 
uniform  electric  field.  They  agree  with  Rabinovitch  and  Zak's  criticism's  of  this 
phenomenon  [Ref.  22],  and  using  Ehrenfest's  theorem,  attempt  to  show  why  Bloch 
oscillations  could  not  occur.  They  consider  the  following  Hamiltonian  of  the  system: 

H  =  p2/2m  +  KV(x)  +  eFx  (2.19) 

where  >^  is  a  constant  which  does  not  affect  the  spatial  distribution  of  the  lattice  potential 
V(x) .  They  attempt  to  apply  Ehrenfest's  theorem  to  (2.19)  [Ref.  23]  and  derive  the 
acceleration  for  the  electron: 

d2(x)/dr  =  (-V[XV(x)  +  eFx])  =  -eF  +  X{-VV)  (2.20) 

By  showing  a  negative  (non-zero)  constant  term  (-eF)  in  the  acceleration,  Churchill  and 

Holmstrom  claim  [Ref.  231  the  velocity  must  decrease  with  time,  and  thus,  Bloch 
oscillations  cannot  occur.  Churchill  and  Holmstrom  did  not  show  the  derivation  of  this 
result.  The  author  of  this  thesis  believes  this  result  may  be  incorrect,  for  they  did  not 
show  that  Ehrenfest's  theorem  was  applied  to  a  wave  function  in  a  Bloch  state  (2.1). 
Also,  they  did  not  specify  the  periodicity  of  the  potential  V(x)  =  V(x  +  a).  Their  result 

appears  so  general  as  to  apply  to  any  particle  accelerated  by  an  electric  field  in  any 
general  potential,  and  it  can  be  shown  there  are  many  instances  where  (2.20)  is  incorrect, 
e.g.,  a  proton  accelerated  by  an  electric  field  under  the  influence  of  gravity. 

In  another  paper  [Ref.  24],  Churchill  and  Holmstrom  build  the  electron  wave 
function  using  eigenstates  assuming,  a  priori,  that  the  energy  states  are  in  the  form  of  a 
Stark  ladder.  They  show  [Ref.  21]  that  their  Bloch-type  wave  functions  do  not  approach 
conventional  zero  field  (F  =  0)  Bloch  states.  Later  analysis  by  Krieger  and  Iafrate  [Ref. 
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23]  show  this  to  be  the  fault  with  their  analysis  and  later  conclusions.  (Note:  Other 
analyses  by  other  researchers  of  Bloch  oscillations  show  that  as  F  -»  0 ,  the  electron 
remains  in  a  Bloch  state  [Refs.  8-11, 14-15,  21-22,  26])  Their  final  conclusion,  again,  is 
that  Bloch  oscillations  do  not  exist,  but  that  "...  accelerated  Bloch  states  can  be  looked 
upon  as  either  a  type  of  standing  wave  or  produced  as  an  interference  of  a  stream  of 
electrons  moving  in  the  +x  direction  and  a  stream  moving  in  the  -x  direction."  [Ref.  24]. 
In  1991  [Ref.  25],  they  attempt  to  refute  Krieger  and  Iafrate's  criticisms  of  previous  work 
[Refs.  23-24],  and  show  that  Bloch  oscillations  exist  in  free  space.  Based  on  their  proof 
[Ref.  25],  they  conclude  that  Bloch  oscillations  cannot  exist  in  a  solid  if  they  exist  in  free 
space.  The  author  of  this  paper  believes  their  conclusion  is  incorrect. 

H.  J.  B.  KRIEGER  AND  G.  J.  IAFRATE 

In  1985,  Krieger  and  Iafrate  [Ref.  26]  summed  up  the  criticisms  [Refs.  14-16,  20- 
25]  of  the  previous  analytical  solutions  to  Bloch  oscillations  by  Zak,  Rabinovitch, 
Churchill  and  Holmstrom.  According  to  Krieger  and  Iafrate  [Ref.  26],  the  criticisms  can 
be  summarized  as  follows: 

1.  The  energy  eigenvalues  of  the  time  independent  Schrodinger  equation  are  not 
quantized  but  are  continuous  with  all  values  of  Er  allowed. 

2.  Since  the  Hamiltonian  is  not  periodic  on  the  boundaries  of  a  (finite)  crystal,  it 
is  not  clear  that  one  can  employ  the  crystal  momentum  representation  of  Houston 
functions  since  Bloch  functions  are  periodic  on  the  boundary,  i.e.,  a  superposition 
of  Bloch  functions  to  represent  the  wave  function  ip  will  automatically  yield  a 

op  which  is  periodic  on  the  boundary,  but  the  solution  of  the  time-dependent 
Schrodinger  equation,  including  the  non-periodic  scalar  potential,  <j)  =  -eEx  may 
not  have  that  property. 

3.  The  crystal  momentum  representation  of  the  operator  x  which  enters  in  the 
calculation  may  not  be  well  defined  because  x$nk  cannot  be  represented  as  a 

linear  combination  of  Bloch  states,  i.e.,  f|x<t>,J2  di  diverges  as  the  crystal 

approaches  the  infinite  extent  in  the  x  direction. 
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Krieger  and  Iafrate  approach  the  solution  of  the  time  evolution  of  electrons  in  an  electric 
field  slightly  different.  Instead  of  using  the  familiar  scalar  potential  for  the  electric  field, 
the  chose  [Ref.  26]  to  represent  it  as  a  vector  potential  which  is  a  function  of  time.  The 
time  dependent  Schrodinger  equation  as  solved  by  Krieger  and  Iafrate  is: 


Hxp(x,t)  = 


(p-(e/c)Ay 


2m 


V(x) 


ty(x,t)  =  ifr 


dq(x,t) 


dt 


(2.21) 


where  A  =  jS(t')dt' ,  and  S{t)  is  an  electric  field  turned  on  at  time  t  =  0 

o 

The  eigenfunctions  [Ref.  26]  are  solved  from  following  equation: 


(p2-(e/c)Ay 


2m 


+  V(x) 


<j>;(*,0  =  s,(04>,'(*>0 


(2.22) 


Using  a  gauge  transformation  [Ref.  26],  the  solution  of  the  eigenfunctions  are: 


ieAxitic, 


^(x,t)  =  eieAX'^n,(x), 


(2.23) 


Expanding  the  wavefunction  ty(x,t)  as  a  sum  of  eigenfunctions  in  Eq  (2.24),  and  solving 
for  the  time  variation  of  the  expansion  coefficients, an(t) ,  of  the  eigenfunctions,  Krieger 
and  Iafrate  [Ref.  26]  reach  the  same  result  as  Houston  [Ref.  7]  did: 


F(t)  _  !r'[E„(*(0)-£„'(*(0)vfc' 

in       r.' 


(2.24) 


Zak  [Ref.  27]  criticized  this  formulation,  and  claimed  this  solution  was  incorrect  due  to 
the  wave  function  repeating  itself  at  every  time  on  every  lattice  cell.  He  concludes  the 
basis  set  (2.24)  is  incorrect,  and  that  Krieger  and  Iafrate  moved  the  problems  of  spatial 
periodicity  into  the  time  domain.  Krieger  and  Iafrate 's  rebuttal  [Ref.  28],  agree  that  their 
basis  states  are  periodic  in  time  at  lattice  sites,  but  disagree  with  his  conclusions 
concerning  this  fact.  They  show  for  any  lattice  site  that  the  wave  function  repeats  itself  at 
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the  Bloch  repeat  time,  not  just  any  time  as  Zak  [Ref.  27]  assumes.  Also,  they  show  in 
Zak's  derivation  that  he  actually  uses  a  different  Hamiltonian  than  they  do  [Ref.  28],  and 
seem  to  overcome  previous  objections  with  boundary  conditions  [Refs.  14-16,  20-25] 
using  their  expansion  on  basis  states  based  on  the  vector  potential  derivation. 

L  L.  ESAKI  AND  R.  TSU 

Esaki  and  Tsu,  in  1970,  predicted  that  Bloch  oscillations  might  be  seen  in  a 
structure  called  a  "semiconductor  superlattice"  [Ref.  3].  Given  the  state  of  art  in 
semiconductor  fabrication,  it  would  be  possible  to  grow  layers  of  different  material  in  a 
way  that  a  periodic  band  structure  would  be  created.  Given  the  probability  of  Zener 
interband  tunneling,  Esaki  and  Tsu  [Ref.  3]  said  that  if  the  scattering  times  were  long, 
Terahertz  (THz)  oscillations  would  be  possible,  and  the  electron  would  undergo 
reflections  at  miniband  boundaries  [Brillouin  Zone].  The  typical  numbers  needed,  as 
estimated  by  Esaki  and  Tsu  for  frequency  vB  =  250  GHz,  F  =  103  V/cm,  and  a  =  10 
nm ,  with  the  scattering  time  greater  than  four  picoseconds.  This  article  was  the  first  in 
which  the  possibility  was  raised  of  creating  a  structure  in  which  radiation  from  an 
electron  undergoing  Bloch  oscillations  might  be  obtained. 

J.  D.  EMLN  AND  C.  F.  HART 

In  1988,  Emin  and  Hart  [Ref.  29]  commented  on  the  time  evolution  of  electrons  in 
a  periodic  potential  under  the  influence  of  an  external  electric  field.  Their  novel  approach 
showed  that  the  applied  electric  field  can  be  broken  up  as  a  sum  of  a  periodic  function 
(sawtooth),  and  non-periodic  (ladder),  and  that  in  their  multiband  solution,  the  ladder 
forms  the  basis  of  the  Stark  ladder  of  multiband  energy  eigenstates  [Ref.  29].  Using  a 
modification  of  Krieger  and  Iafrate's  approach  [Ref.  26],  they  derive  an  equation  of 
motion  where  their  Bloch  eigenstates  are  a  function  of  the  periodic  portion  of  the  electric 
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field  [Ref.  29].    Since  then,  the  Emin-Hart  approach  has  been  shown  by  several  authors 
to  be  incorrect  [Refs.  30-32]. 

K.  M.  LUBAN  AND  A.  BOUCHARD 

In  1993,  Luban  and  Bouchard  [Ref.  5],  using  highly  accurate  numerical 
techniques,  simulated  the  motion  of  an  electron  wavepacket  in  a  semiconductor 
superlattice  under  the  influence  of  an  electric  field.  The  advantage  of  their  approach  was 
the  numerical  simulation  contained  all  powers  of  the  Hamiltonian  containing  all  energy 
bands,  instead  of  single  band  Hamiltonian  approximations  used  by  previous  analytical 
approaches.  Coinciding  with  recent  four  wavemixing  experiments  which  detected  Bloch 
radiation  for  the  first  time  from  superlattices  [Ref.  4],  Luban  and  Bouchard's  numerical 
simulations  clearly  demonstrated  (absent  scattering  events)  that  an  electron  will  undergo 
Bloch  oscillations  in  a  semiconductor  superlattice,  and  for  the  first  time,  showed  that 
Bloch  oscillations  are  a  bona-fide  aspect  of  the  dynamics  of  independent  electrons  in  a 
solid. 
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m.  NUMERICAL  METHODS 

A.  TIME  DEPENDENT  SCHRODINGER  EQUATION  (TDSE) 

Solving  the  Time  Dependent  Schrodinger  Equation  is  essential  to  the  modeling  of 
Bloch  oscillations.  The  general  one  dimensional  Schrodinger  equation  is: 

+  V(x,t)y(x,t)  =  ih — - ,  (3.1) 


2/n       dx2  dt 

where  ty(x,t)  is  the  wave  function  of  the  electron,  and  V(x,t)is  the  potential  energy 

function  of  the  electron.  We  will  consider  a  system  in  which  V(x,  t)  is  independent  of 

time,  thus  making  V  a  function  of  position  only.  The  Hamiltonian  of  a  system  is  defined 

as  H  =  T  +  V ,  where  T  is  the  kinetic  energy  operator,  and  V  is  the  potential  seen  by  the 

electron.  The  Schrodinger  equation  can  be  cast  in  terms  of  the  Hamiltonian  operator.  If 

h  d 

the  momentum  operator  p  is  given  by ,  then  the  Schrodinger  equation  can  be  cast 

i  dx 

2 

in  the  form  of  an  operator  equation,  with  H  =  ^—  +  V .  The  time  dependent  Schrodinger 

2m 

equation  reads: 

Hxp(x,t)  =  ih-^-  (3.2) 

ot 

To  solve  this  numerically,  one  must  discretize  (3.2).  We  will  discuss  the  discretization 
process  in  the  following  sub-sections. 
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1.  Discretization  of  the  Wave  Function 

First,  we  represent  the  continous  variable  x  as  a  discrete  variable,  x  =  x0+  nAx , 
where  Ax  is  the  mesh  size,  x0  is  the  starting  value,  and  n  is  any  integer.  Also,  we 
represent  time  as  a  discrete  variable,  t  =  tQ  +  mAt ,  where  At  is  the  time  mesh  size,  t0  is 
the  initial  time,  and  m  is  any  integer.  The  initial  points,  x0,t0  are  both  set  equal  to  zero 
to  simplify  the  scheme.  Also,  this  form  of  discretization  for  time  allows  the  efficient  use 
of  a  loop  in  a  computer  program.  The  form  of  our  sampled  wave  function  is  a  column 
vector,  with  each  sample  being  taken  at  a  mesh  point,  such  that  for  any  time  t ,  the  wave 
function  will  look  like: 


Vm*(x)  = 


ty(Ax) 
U;(2Ax) 

\\)((N  -  l)Ax) 

\\>(NAx) 


(3.3) 


where  N  is  the  total  number  of  mesh  points  in  the  discretization.  Using  a  loop  over  time, 
the  computer  program  generates  a  new  op  vector  for  each  time  step. 

2.  Discretization  of  the  Hamiltonian 


*2       -2 

n      d 

The  Hamiltonian  has  two  terms:  the  kinetic  energv  operator  T  = ,  and 

2m  dx~ 

the  potential  V(x) .  We  must  find  a  suitable  discrete  representation  for  each  term.  The 
discrete  representation  of  the  potential  is  easy. 
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Being  a  function  of  x ,  it  takes  on  the  same  form  as  the  wave  function,  i.e.,  a  column 
vector  such  that: 


V(*)- 


V(Ax) 
V(2Ax) 

V((N-l)Ax) 

V(NAx) 


(3.4) 


The  discretization  of  the  kinetic  energy  operator  is  more  involved.  A  suitable  discrete 
representation  for  the  second  spatial  derivative  must  be  found.  To  approximate  it,  we 
look  at  the  Taylor  series  expansion  for  arbitrary  functions  f(x  +  Ax)  and  f(x  -  Ax) : 


1  df(x)  1  dzf(x)      , 

f(x  +  Ax)=  f(x)  +  -  -L^J-Ax  +  -     -^Ar^Ax2- 


1!    dx 


2!    dx' 


(3.5) 


and 


1  df(x)  1  d2f(x)      , 

f(x-Ax)=  fix) L^J-Ax-¥ J-ArLAx2+. 

M  '     JK  '    1!    dx  2!    dx2 


(3.6) 


Ignoring  terms  higher  than  order  two,  and  taking  the  sum  of  (3.5)  and  (3.6),  we  find  that 


d2f(x)      f(x  +  Ax)  +  f(x  -  Ax)  -  2f(x) 


dx' 


Ax' 


(3.7) 


Letting  the  total  derivative  go  to  a  partial  derivative  in  x ,  and  letting  ip  be  our  arbitrary 
function  / ,  the  kinetic  energy  operator  T  acting  on  ip  for  an  arbitrary  mesh  point 
nAxis: 


2       2                                 2 
jip(/iAx)  = ^-[^((/i+l)Ax)  +  T|j((«-l)Ax)-2ip(nAr)].       (3.8) 
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Therefore,  the  Hamiltonian  acting  on  xp  at  an  arbitrary  mesh  point  n Ax  has  the  form: 


AAj)(/iAx)  = 


2m  Ax 


-[tj>((n  +  l)Ax)  +  \J>((/i  -  l)Ax)  -  2\J>(rtAx)]  +  v(«Ax)tp(nAx)  .  (3.9) 


If  our  mesh  runs  from  n  =  0,1,2...  Af ,  in  solving  the  the  Hamiltonian,  it  arises  that  values 
for  \j;(-lAx)  and  ty((N  +  l)Ax)  are  needed.  These  exist  outside  the  mesh  system  we  are 

looking  at,  and  therefore  are  simply  set  to  zero.  This  is  tantamount  to  assuming  the 
potential  suddenly  becomes  infinite  outside  the  mesh  system.  This,  as  a  practical 
consequence,  forces  one  to  construct  the  spatial  system  sufficiently  large  so  that,  over  the 
time  scales  of  interest,  the  wave  functions  do  not  reflect  from  these  "hard  walls." 

The  effect  of  the  discretization  is  thus  to  replace  a  single  differential  equation  with 
a  set  of  N  coupled  equations.  The  form  of  the  Hamiltonian  is  an  N  X  N  matrix: 


//  = 


2a  +  V 

a 

0 

... 

0 

a 

-2a  +  V 

a 

0 

: 

0 

*  ^ 

'-. 

■. 

0 

: 

0 

a 

-2a  +  V 

a 

0 

0 

a 

-2a  +V 

(3.10) 


where  a  =  - 


fr 


2m  Ax" 


.  It  is  seen  that  the  Hamiltonian  is  a  tridiagonal  matrix,  with  only 


the  main  diagonal  and  the  first  upper  and  lower  diagonals  having  non  zero  terms.    To 
reduce  the  amount  of  space  required  for  computer  memory,  the  computer  program  in 
Appendix  A  uses  a  sparse  representation  for  storing  the  Hamiltonian.  Instead  of  an 
NX  N  matrix,  we  use  an  N X  3 matrix  storing  only  the  non-zero  terms: 


H  = 


0  -2a  +  V  a 

a  -2a  +  V  a 

a  -2a  +  V  a 

a  -2a  +  V  0 


(3.H) 
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B.  SOLUTIONS  OF  THE  TIME  DEPENDENT  SCHRODINGER 
EQUATION 

Equation  (3.2)  is  a  first  order  differential  equation  of  the  time  variable.  The 
solution  of  this  first  order  differential  equation  has  a  unique  solution,  given  in  the 
Schrodinger  representation  by: 

iHi 

q(x,t)  =  e~ n  y(x,0)  (3.12) 

To  perform  practical  computations  of  ty(x,t) ,  we  must  know  two  things.  First,  the  initial 
state  of  the  wave  function  must  be  known  or  guessed  for  time  t  =  0 .  Secondly,  an 

approximation  for  e     h   must  be  made.  This  value  is  not  only  complex,  but  it  is  also 
strictly  unitary,  and  any  approximation  must  maintain  unitarity  for  it  to  be  considered 
valid.    Unitarity  means  that  the  normalization  of  the  wave  function  is  constant  at  all 
times,  and  the  wave  function  for  any  time  t  must  satisfy  the  following  integral  equation: 

00 

f\b'(x,t)ty(x,t)dx  =  l  (3.13) 

—x 

In  the  following  sub-sections,  we  will  discuss  approximation  methods  and  their  relative 
strengths  and  weaknesses. 

1.  Power  Series  Expansion 

Numerically  calculating  the  time  evolution  operator  e  h  is  not  easy.  One  could 
use  a  power  series  expansion  to  represent  the  operator.  For  an  arbitrary  operator  A ,  the 
power  series  expansion  of  eA  is: 
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A2      A3 
eA  =l  +  A  +  —  +  —  +...  (3.14) 

2!      3! 

Based  on  this  result,  the  power  series  expansion  of  the  time  evolution  operator  is: 

-m          iHt     H2t2     iHY  „_ 

e   n   =1 =-  + r-+...  (3.15) 

n     2n2     en3  v     ' 

A  truncation  of  the  expansion  can  be  implemented  easily  numerically,  since  it  is  only  a 
finite  sum  of  products  of  terms.  The  restriction  that  any  approximation  maintain  unitarity 
and  satisfy  (3.13)  prevents  the  use  of  it  in  program.  Any  truncation  of  this  infinite  series 
removes  the  property  of  unitarity  from  it,  and  since  a  computer  cannot  generate  an 
infinite  amount  of  terms,  any  solution  generated  from  the  truncation  would  be  invalid. 

2.  Magnus  Expansion 

In  1954,  Wilhem  Magnus  derived  an  expansion  which  can  be  used  to  represent  the 
time  evolution  operator.  Using  the  Baker-Hausdorff  theorem,  Magnus  derived  a  solution 
to  the  first  order  differential  equation  [Ref.  33]: 

4-Y(t)  =  A(t)Y(t)  (3.16) 

at 

where  Y(t)  is  of  the  form: 

Y(t)  =  eQ(,)Y(0)  (3.17) 

and  Q(t)  is: 


Q(t)  =fA(xyh  +  \  f A(T),}A(o)do]dx+. . .  (3.18) 

o  2  "o  c  J 

where  the  brackets  indicate  commutators.  This  result  was  independently  derived  by 
Pechukas  and  Light  [Ref  34].  Applying  the  Magnus  expansion  requires  a  very 
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sophisticated  numerical  integration  routine,  and  requires  a  great  deal  of  error  analysis  and 
monitoring  of  the  results  to  ensure  numerical  stability  and  accuracy  [Ref.  35].  Because  of 
its  complexity,  this  method  was  not  chosen. 

3.  Explicit  Alogrithm  for  Solving  the  TDSE 

In  1991,  P.  B.  Visscher  developed  an  explicit  numerical  algorithm  to  solve  the 
time  depedent  Schrodinger  equation  [Ref.  36].  Starting  with  equation  (3.1),  Visscher 
rewrote  it  in  terms  of  solutions  of  the  real  and  imaginary  portions  of  the  wave  function 
such  that: 


dty(x,t) 
VrV      '  =  Hty.(x,t)  (3.19) 

at  ' 

and 


*£2~tt,r(*0  (3.20) 

at 

where  tyr(x,t)  is  the  real  portion  of  the  wave  function,  and  tp.(;t,f)is  the  imaginary 

portion  of  the  wave  function  [Ref.  36].  He  claims  this  approach  is  the  same  as  calculating 
a  trajectory,  and  thus  in  the  discretization  of  each  equation,  each  portion  is  defined  at 
staggered  times.  For  the  real  portion  of  the  wave  function  tyr(x,  t) ,  it  is  defined  at  times 

equal  to  t  =  0,  At, 2  At,. . . ,  while  the  imaginary  portion  is  defined  at  times 

t  =  05At,15At,25At,...  [Ref.  36].  The  discretization  of  equations  (3.19)  and  (3.20)  is: 

yr(x,t  +  05 At)  =  \pr(x>t-  05At)  +  AiH\\>i(x,t)  (3.21) 


anc 


qt(x,t  +  05At)  =  u;:(;t,f  -  05At)  -  AtHyr(x,t)  (3.22) 
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To  start  the  iteration,  initial  values  for  tyr(x,0)  and  ii),(;c,0.5Ar)  must  be  obtained,  and 

time  must  start  at  t  =  05 At  [Ref.  36].  Visscher  also  points  out  that  because  the  real  and 
imaginary  portions  of  the  wave  function  exist  at  different  time  steps,  calculating  the 
overall  probability  density  could  be  problematic.  A  set  of  equations  which  conserve 
probability  for  integer  and  half  integer  units  of  time  are: 

\ty(x,tf  =  op  r(x,tf  +  ip  t{x,t  +  05At)\p  ,{x,t  -  05  At)  (3.23) 

for  integer  time  intervals,  and  for  half  integer  time  intervals: 

\\p(x,tf  =  y  r(x,t  +  05At)ty  r(x,t  -  05At)  +  ip  i(x,t)2  (3.24) 

This  scheme  is  second  order  accurate,  and  conserves  probability,  satisfying  equation 
(3.13).  The  claim  of  Visscher  is  this  scheme  is  three  times  faster  than  the  actual  scheme 
used  by  the  author  of  this  thesis  (Cayley  method).  This  paper  was  discovered  in  the 
process  of  the  thesis  write-up,  and  therefore  this  method  and  its  author's  claims  could  not 
be  compared  to  the  Cayley  method.  If  Visscher's  claim  is  true,  then  it  would  be  an 
outstanding  method  to  use,  and  would  be  an  excellent  model  for  future  students  to  use. 

4.  Cayley  Hamilton  Approximation 

The  approximation  actually  used  for  the  time  evolution  operator  is  the  Cayley- 
Hamilton  approximation  [Refs.  37,38].  It  is  Hermitean,  and  maintains  unitarity,  as  well 
as  has  good  properties  for  keeping  error  low.  The  Cayley  approximation  of  the  time 
evolution  operator  is  given  by: 


m      1-^iHAt 
e    *    . 28 (3.25) 
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which  has  error  on  the  order  of  (At)  .  This  error  can  be  seen  when  compared  to  the 


iHM 


power  series  expansion  of  e    h   .  Taking  the  first  four  terms  in  the  power  series 
expansion  of  this  time  evolution  operator  from  (3.14): 


-—■           iHAt     (HAt)2       (HAt)3 
e    *    =1 ^ j-  +  i± /-+...,  (3.26) 

«        2/r  6/zJ 

while  the  expansion  for  the  approximation - is: 

1  +  jt  iHAt 

In 


l-^iHAt  _l     iHAt     {HAt)2  ^i(HAt)3  ; 
1  +  ^iHAt  A  2h2  A%3 

these  two  expansions  differ  beginning  only  in  the  third  term,  so  the  error  will  be  third 
order.  The  Cayley  approximation  to  the  original  solution  is: 


1- %  iHAt 

1  +  ^-iHAi 


ty(x,t  +  At)  =  ; — i  — :  ty(x,  t)  (3.28) 

Finally,  the  computer  will  iterate  the  following  algorithm: 

(1  +  jKiHAt)^(x,t  +  At)  =  (1  -J£iHAt)y(x,t)  (3.29) 

In  numerical  analysis,  this  type  of  solution  is  called  an  "implicit"  scheme  or  a  Crank- 
Nicholson  scheme  [Ref.  38]. 

C.  SOLVING  THE  CAYLEY  APPROXIMATION 

Equation  (3.29)  has  the  form: 

M\{x,  t  +  At)=  My(x,  t)  (3.30) 


25 


where  M  =  I  +  —  H ,  and  /  is  the  identity  matrix,  and  H  is  an  NXN  matrix. 

2% 

Solving  for  the  wave  function  at  the  next  time  step  is  the  ultimate  goal  in  order  to  model 
the  time  varying  motion  of  the  electron.  There  are  several  methods  of  solving  equation 
(3.29)  available,  and  the  methods  will  be  discussed  in  the  following  sub-sections. 

1.  Solution  by  Inversion 

The  most  direct  solution  of  equation  (3.30)  for  ty(x,t  +  At)  is: 

y(x,t  +  At)  =  (M')~i  My(x,t)  (3.31) 


-i 


but  it  is  numerically  intensive.   M   is  a  tridiagonal  matrix,  while  the  inverse,  (M  )    ,  is 

a  full  matrix,  and  must  be  solved  before  proceeding  to  solve  equation  (3.31).  The  inverse 
can  be  obtained  by  several  methods  (Gaussian  reduction,  Gauss  Jordan  elimination, 
etc.),  but  the  ultimate  solution  of  equation  (3.31)  takes  an  order  of  N3  operations  to 

solve  [Ref.  38]  if  (M*)     and  M are  N  X  TV  matrices.  There  is  a  more  elegant  method 

which  takes  advantage  of  the  tridiagonal  nature  of  M  and  M  * ,  and  reduces  the  number 
of  operations  to  the  order  of  N .  This  is  the  L-U  decomposition  method,  and  will  be 
discussed  in  detail  next. 

2.  L-U  Decomposition 

This  method  of  solving  a  linear  equation  was  taken  from  Numerical  Recipes  in  C 
[Ref.  ].  Equation  (3.30)  has  the  form  A  •  x  =  b ,  where  b  is  an  TV  X  1  column  vector  and 
is  the  solution  of  [/  -  -^  H]ip(x,  t).  A  =  [/  +  4j£  H]  =  M* ,  and  x  is  what  we  are  solving 

for:  \l)(x,t  +  At). 
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In  order  to  solve  this  by  L-U  decomposition,  we  break  A  up  into  two  matrices  such  that: 


A-x  =  (L-U)-x  =  L-(U-x)  =  b. 


(3.32) 


We  solve  this  by  back  substitution,  first  solving  for  an  intermediate  vector  v ,  such  that 
L-  v  =  b,  and  then  back  solve  U  •  x  =  v .  This  involves  solving  a  set  of  TV  linearly 
independent  equations.  In  decomposing  a  tridiagonal  matrix,  L   takes  on  the  form: 


L  = 


an 

0 

... 

o  ■ 

a21 

a22 

I 

0 

"• 

'• 

: 

aAf-2,Af-l 

aAM,.V-l 

0 

0 

... 

0 

V-N-ljV 

aNN 

(3.33) 


where  a^  are  the  coefficients  of  the  L  matrix.  In  order  to  conserve  computer  storage 

memeory,  the  program  in  the  Appendix  stored  the  L  matrix  in  an  N  X  2 ,  which  only 
held  the  non-zero  elements,  so  that  L  is: 


L  = 


a, 


a 
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0 
(X-,-, 


CtAT-2^-l       aAf-UV-l 


(3.34) 


Likewise,  the  U  matrix  was  handled  in  a  similar  fashion.  The  form  of  the  U  matrix  when 
derived  from  a  tridiagonal  matrix  is: 


U  = 


■p,l 

Pl2 

0 

... 

0 

0 

P22 

\% 

•. 

0 

Pa 

-!JV-: 

Pj\'-l,N 

0 

0 

Paw 

(3.35) 


where  8,  are  the  coefficients  of  the  U  matrix. 

27 


Again,  to  conserve  space,  U  was  actually  stored  in  an  N  X  2  matrix  and  takes  almost  the 
same  form  as  L  : 


U  = 


fcl 

fta 

022 

^23 

Pn-ijj-i 

Pn-im 

0 

h 

(3.36) 


The  method  of  solution  is  straight  forward.  First,  all  the  coefficients  of  the  L  and 

U  matrices  must  be  solved.  Then  equation  (3.32)  can  be  solved  in  two  steps.  The  first 

step  is  to  solve  for  L  •  y  =  b ,  which  can  be  done  as  follows: 


tt 


A 
a,, 


y,  = 


ct. 


for  i  =  2,3,. ..,JV 


(3.37) 


This  procedure  is  accomplished  very  quickly  due  to  the  tridiagonal  nature  of  the  L 
matrix.  Then  to  solve  for  x ,  the  final  answer,  it  is  solved  by  back  substitution  as  follows: 


jV  — 


>V 


'AW 


1 

x,  =  — 

P4 


*- 2  to 


for  /  =  7V-l,iV-2,...,l 


(3.38) 


Once  this  answer  for  the  wave  function  for  the  new  time  step  is  obtained,  it  is  multiplied 
by  M  to  form  a  new  b  vector,  and  the  process  is  repeated  to  solve  for  the  wave  function 
at  the  next  time  step.  This  process  is  repeated  in  a  loop  over  the  time  of  the  simulation. 
Results  of  the  program  will  be  discussed  in  the  next  chapter. 
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IV.  RESULTS  AND  CONCLUSIONS 

The  computer  code  in  the  Appendix  has  been  under  development  for  nearly  a 
year.    The  solution  of  the  time  dependent  Schrodinger  equation  using  the  Cayley 
approximation  sounds  easy  to  accomplish  numerically,  but  has  been  challenging  to 
successfully  accomplish.  C  language  was  chosen  due  to  the  local  support  and  our 
familiarity  with  it.  Since  the  electron  wave  function  is  complex,  solving  the  time 
dependent  Schrodinger  equation  numerically  required  the  use  of  a  user  defined  structure. 
All  variables  were  defined  of  type  complex,  with  a  real  and  imaginary  part  of  double 
precision.  The  first  efforts  in  developing  the  code  were  writing  reliable  functions  to  add, 
subtract,  multiply,  divide,  and  conjugate  complex  numbers.  Tnese  tasks  were  verified 
using  simple  algebraic  equations  that  the  author  programmed  and  knew  the  results  ahead 
of  time.  The  functions  performed  as  expected. 

The  next  step  was  to  write  a  routine  that  would  perform  an  L-U  decomposition  on 
an  N  X  N  matrix.    The  concepts  were  well  laid  out  in  Numerical  Recipes  in  C  [Ref.  38], 
but  tailoring  it  a  tridiagonal  matrix  was  more  difficult  than  expected.  After  several 
attempts,  C  code  was  written  to  perform  the  L-U  decomposition,  and  results  of  our 
decomposition  routine  on  a  general  NXN  matrix  were  verified  against  an  identical 
matrrix  using  MATLAB's  L-U  decomposition  function.  The  decomposed  matrices  were 
identical.  Finally,  the  back  substitution  routine  was  written  to  solve  the  time  dependent 
Schrodinger  equation  (3.29)  using  the  Cayley  approximation.  To  ensure  our  code  was 
accurate,  our  first  model  was  a  one  dimensional  free  particle,  with  V  =  0 .  Tne  wave 
packet  for  this  model,  and  all  subsequent  simulations  was  a  normalized  Gaussian  wave 

packet  yp(x)  =  Ce  2a  .  where  C  is  a  normalization  constant,  and  a  controls  the  width  of 

the  packet.  The  free  particle  model  behaved  as  expected,  with  the  wave  packet  spreading 
out  in  both  the  +x  and  -x  directions.    To  ensure  our  numerical  model  was  correct, 
another  system  simulated  was  an  electron  wave  packet  in  a  harmonic  oscillator  potential 
V  =  \  moi2x2 .  The  response  of  the  particle  should  be  to  diffuse  out,  and  gather  itself  back 
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at  exact  time  intervals  of  2ji/o)  .  The  model  again  showed  the  wave  packet  evolve  as 
expected,  and  repeat  at  exact  time  intervals  of  2ji/o)  . 

Confident  that  the  model  to  numerically  simulate  the  time  dependent  Schrodinger 
equation  was  accurate,  electron  wave  packets  in  periodic  potentials  with  an  external 
electric  field  applied  were  modeled.  The  simulations  were  ran  either  on  a  Sun 
Sparcstation  20,  or  a  Cray  C-90  super-computer  at  Wright  Patterson  AFB.  For  problems 
modeling  Bloch  oscillations,  the  Sparcstation  typically  took  12  hours  to  run,  while  the 
super-computer  would  take  eight  hours.  A  variety  of  results  were  generated,  and  a 
representative  sample  is  presented  and  discussed  below: 

A.  NON  LOCALIZED  WAVE  PACKET  IN  GENERAL  POTENTIAL 

Figure  (1)  is  a  picture  of  an  electron  clearly  undergoing  Bloch  oscillations. 


Figure  (1).  Electron  Undergoing  Bloch  Oscillations 

For  all  simulation  output,  the  vertical  axis  represents  time,  with  the  electron  wave  packet 
starting  at  the  top  of  the  figure  working  its  way  down  in  time,  and  the  horizontal  axis  of 
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ihe  figure  is  space.  The  color  bar  at  the  bottom  of  the  figure  shows  the  relative 
probability,  with  zero  probability  being  green,  and  maximum  probability  being  yellow, 
with  a  continuous  color  scale  showing  probabilities  in  between.    The  horizontal  scale  is 
10  nm,  and  the  amplitude  of  oscillation  being  nearly  that  large.  For  the  system  in  Figure 
(1),  the  potential  function  is  given  by  V  =  0^|cos(2x)j ,  with  a  =  jt/2  nm  being  the  repeat 

distance.  The  applied  electric  field  strength  is  F  =  105  V/cm .  Given  the  previous 

parameters,  the  Bloch  repeat  time  is  tb  =  -^  =  264  fs .  The  code  was  run  for  16000  time 

steps,  v/hich  generated  three  oscillations.  It  is  clear  the  center  of  mass  of  the  electron  is 
oscillating  back  and  forth,  and  thus,  we  would  expect  to  radiate  like  a  classical  dipole. 
No  other  dynamic  electron  effects  are  seen  here,  so,  by  accident,  we  chose  a  system  that 
would  undergo  perfect  Bloch  oscillations  instead  of  executing  interband  transitions  or 
remaining  in  a  stationary  state.  Without  calculating  the  energy  bands  and  eigenvectors 
and  composing  the  inital  state  of  the  electron  as  a  superposition  of  eigenfunctions  of  a 
single  band  (Wannier  functions),  the  behavior  could  not  be  exactly  predicted  prior  io 
running  the  code. 
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B.  LOCALIZED  ELECTRON  IN  A  GENERAL  POTENTIAL 


Figure  (2)  is  a  picture  of  an  electron  undergoing  a  "breathing  mode"  of  oscillation. 


Figure  (2).  Electron  in  a  "Breathing  Mode" 

The  system  is  exactly  the  same  as  previous,  with  the  only  difference  being  that  the 
electron  wave  packet  is  localized  in  a  well.  We  did  this  by  reducing  o ,  which  controls 
the  width  of  the  wave  packet.  Instead  of  executing  center  of  mass  oscillations,  the 
electron  wave  packet  "breathed"  in  and  out  at  the  Bloch  frequency.  The  electron 
preferentially  breathed  in  the  direction  of  lower  energy.  The  code  was  run  for  16000 
time  steps,  showing  three  breathing  modes.  This  phenomenon  would  not  radiate, 
because,  as  seen  in  the  picture,  the  average  center  of  mass  does  not  appreciably  move 
over  time.  This  phenomenon  was  actually  predicted  by  Zener  in  1933  [Ref.  2],  and  his 
explanation  is  the  electron  is  localized  near  the  Brillouin  zone  boundary,  and  it  only 
transitions  between  k  =  +  f  and  k  -  -f  in  k  space  by  Bragg  reflection,  thus  executing  a 

"breathing"  in  real  space. 
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C.  ELECTRON  REFLECTIONS  AT  BOUNDARIES 

Figure  (3)  is  an  example  of  electron  wave  packet  reflection  off  the  boundaries  of 
the  system  investigated: 


Figure  (3).  Electron  Reflection  of  Boundaries 

In  the  process  of  discretizing  the  time  dependent  Schrodinger  equation,  it  imposes  finite 
boundaries  on  the  system.  The  effect  of  reaching  a  boundary,  as  seen  in  Figure  (3),  is 
quite  dramatic.  The  wave  reflects  perfectly  off  the  boundary,  and  begins  to  interfere  with 
the  rest  of  the  wave  function  executing  a  Bloch  "breathing"  oscillation.  After  several 
reflections,  the  "breathing"  becomes  completely  distorted,  as  the  wave  function  folds 
over  on  the  right  side.  For  future  simulations,  a  test  could  be  developed  to  determine  if 
the  wave  function  was  near  a  boundary,  and  the  test  would  stop  the  simulation  if  a 
reflection  occurred. 
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D.  KRONNIG-PENNY  SIMULATIONS 

The  Kronnig-Penny  one  dimensional  square  wave  potential  is  an  excellent  model 
to  simulate  the  conduction  band  edge  of  a  general  semiconductor  superlattice.  The 
superlattice  simulated  is  a  GaAs/Al  0  3  Ga  0  7  As  superlattice,  with  barrier  height  equal  to 
0.243  eV .  In  the  following  sections,  we  will  discuss  the  effect  of  varying  different 
parameters,  and  their  effect  of  the  motion  of  the  electron  in  the  lattice. 

1.  Non-localized  Electron  in  a  Superlattice 

Figure  (4)  is  a  picture  of  an  electron  in  a  superlattice  potential: 


Figure  (4).  Non-localized  Electron  in  a  Superlattice 


The  potential  function  is  given  by  V  =  0.243  eV  for  the  Al  0  3  Ga  0  7  As  barriers,  V  =  0 
eV  for  the  GaAs  barriers,  and  a  =  10  nm  being  the  repeat  distance.  The  AlGaAs  barriers 
were  5.1  nm  thick,  and  the  GaAs  barriers  were  4.9  nm  thick.  The  applied  electric  field 
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strength  is  F  =  104  V/cm .  Given  the  previous  parameters,  the  Bloch  repeat  time  is 
xb  =lk  =  827/5,  and  the  amplitude  of  oscillation  was  approximately  20  nm.  The 
simulation  was  run  for  17000  time  steps,  giving  two  oscillation  periods  to  observe.  The 
electron  wave  packet  was  approximately  250  angstroms  wide.  It  is  seen  part  of  the 
electron  is  "ripped  away"  and  undergoes  Bloch  oscillations,  and  parts  of  it  bifurcate  and 
tunnel  through  the  barriers,  remaining  relatively  stationary.  This  electron  is  clearly  in  a 
superposition  of  a  dynamic  and  stationary  state.  The  electron  is  undergoing  Bloch 
oscillations,  but  it  could  also  be  transistioning  energy  bands  or  remaining  in  the  ground 
state.  Without  calculating  energy  band  transition  probabilities,  it  is  impossible  to  tell  in 
this  calculation.  That  could  and  should  be  built  into  a  more  sophisticated  model. 

2.  Effect  of  Varying  Repeat  Distance 

Figure  (5)  is  a  picture  of  an  electron  in  a  superlattice  with  half  the  repeat  distance 
of  the  electron  in  section  1: 


Figure  (5).  Electron  in  Superlattice  (Differing  Repeat  Distance) 
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The  repeat  distance  was  changed  from  the  system  in  section  1  from  10  nm  to  5  nm. 
Keeping  the  electric  field  constant,  this  changes  the  repeat  time  to  xB  =  -^-  =  4135  fs . 

The  simulation  was  run  for  17000  time  steps,  showing  four  oscillation  periods.  The 
result  is  much  the  same  as  the  previous  two  before,  except  the  electron  tunnels  through 
more  barriers,  and  less  of  it  is  available  for  Bloch  oscillations,  which  are  faintly  present. 

3.  Effect  of  Varying  Barrier  Thickness 

Figure  (6)  shows  the  effect  of  changing  the  system  in  section  2  by  reducing  the 
width  of  the  AlGaAs  barriers,  and  increasing  the  GaAs  well  width  such  that  the  5  nm 
repeat  distance  is  maintained: 


Figure  (6).  Electron  in  a  Superlattice  (Varying  Barrier  Width) 


The  electron  completely  tunnels  through  the  barriers,  and  it  appears  the  electron  is  not 
undergoing  Bloch  oscillations  at  all.  This  may  be  because  the  barriers  are  relatively  thin 
(1.5  nm),  and  repeat  very  closely  to  one  another  (5  nm). 
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4.  Effect  of  Varying  Electron  Mass 

Figure  (7)  is  a  picture  of  the  effect  of  changing  the  electron  mass  of  the  system 
described  in  section  3: 


Figure  (7).  Electron  in  a  Superlattice  (Varying  Electron  Mass) 

The  mass  of  the  electron  in  the  system  in  section  5  was  changed  to  0.85  times  the  bare 
electron  mass.  The  effect  of  making  the  electron  "lighter"  was  to  make  Bloch  oscillations 
somewhat  apparent,  where  as  they  were  not  in  the  previous  simulation.  This  represents  a 
physical  effect.  In  a  semiconductor,  the  interaction  of  an  electron  with  each  material  can 
be  modeled  by  changing  the  electron  mass  to  an  effective  mass.  While  the  actual 
effective  masses  of  GaAs  and  AlGaAs  were  not  used,  the  simulation  is  effective  in 
demonstrating  the  effect  of  the  effective  mass  approximation  in  semiconductors,  and 
shows  that  semiconductors  with  effectively  "light"  electrons  are  ideal  materials  to 
stimulate  Bloch  oscillations  where  as  they  might  not  be  present  on  other  materials. 


37 


E.  SUMMARY 

We  have  demonstrated  that  Bloch  oscillations  are  a  genuine  phenomenon  of  an 
electron  in  a  periodic  potential.  The  rich  variety  of  dynamic  phenoma  of  Bloch  electrons 
is  a  great  source  for  further  research.  There  are  a  wide  variety  of  potential  applications 
for  stable  GHz  and  THz  oscillators  for  high  speed  integrated  circuits  and  Rf  sensors. 
More  sophisticated  models  of  this  effect  can  be  made  such  as  including  the  effective  mass 
variation  of  the  layers,  calculating  the  energy  eigenvalues  and  vectors  to  start  the  initial 
state  in  a  single  energy  band,  calculating  the  effect  of  the  electron  motion  on  the  potential 
using  the  Poisson  equation,  and  using  that  potential  in  the  Schrodinger  equation  (self- 
consistent  Poisson-Schrodinger  equation),  incorporate  scattering  mechanisms  such  as 
electron-electron  scattering,  electron-hole  elimination,  electron  phonon/photon  scattering, 
and  a  variety  of  other  mechanisms  present  in  the  solid.  Including  these  effects  adds  to  the 
complexity  of  the  model  and  code,  and  would  require  the  immense  resources  of  a  super 
computer,  but  accurately  modeling  this  effect  could  reveal  new  and  useful  applications 
for  future  devices. 
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APPENDIX:  CAYLEY  METHOD  IN  ANSI  C 

#include  <stdio.h> 
#include  <stdlib.h> 
#include  <math.h> 
#definePI  3.14159265 

struct  complex  { 
double  r; 
double  c; 

}; 

struct  complex  Cadd(struct  complex  x,  struct  complex  y); 
struct  complex  Csub(struct  complex  x,  struct  complex  y); 
struct  complex  Cmult(struct  complex  x,  struct  complex  y); 
struct  complex  Cdiv(struct  complex  x,  struct  complex  y); 
struct  complex  Ceq(struct  complex  x); 
struct  complex  Cconj  (struct  complex  x); 


struct  complex  CI,  C2,  M[10002][4],  Ms[10002][4],  psi[10002]; 

struct  complex  npsi[  10002],  magpsi[10002],  one,  eye,  sum,  temp; 

struct  complex  L[10002][3],  U[10002][3],  y [10002],  conj[  10002]; 

struct  complex  zero,  tempi,  V[10002],  b[  10002],  psit[10002],  magipsi[10002]; 

struct  complex  nconj[  10002],  nmagpsif  10002]; 

double  dt,dx,  hbar,m,cl,  c2,  tconst,  sigma,  e,  F; 
int  i,j,N; 
long  loop; 

mainQ 

/*  Time  Dependent  Schroedinger  Equation  -  using  Cayley  method  */ 

/*  L-U  Decomposition  program,  code  based  on  procedure  in  Numerical  Recipes*/ 

/*  in  C,  pp  43-44  */ 

/*  L  is  alpha  units  */ 
/*  U  is  beta  units  */ 
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{ 

N=10001; 

m=9.1e-31;  /*  Mass  of  electron  in  kg  */ 

dt=0.05;  /*  Time  step  in  femto  seconds  */ 

dx=0.01; 

hbar=0.6582;    /*  Plancks  constant  divided  by  2  pi  in  eV  fs*/ 

tconst=0.03810;     /*  hbar  squared  divided  by  twice  mass  of  electron  eV  nm2  */ 

c2=tconst/dx/dx;   /*  Constant  that  wraps  up  all  constants*/ 

cl=0.5*dt/hbar;     /*  second  constant  that  wraps  up  all  constants  */ 

sigma=4000.0; 

e=1.0;        /*  Charge  of  electron  */ 

F=le-2;  /*  Field  strength  in  Volts/nm  */ 

Cl.r=0.0; 
Cl.c=1.0*cl; 

C2.r=1.0*c2; 

C2.c=0.0; 

zero.r=0.0; 
zero.c=0.0; 


one.r=1.0;  /*  One  -  real  number  */ 

one.c=0.0; 

eye.r=0.0;  /*  Imaginary  i*/ 

eye.c=1.0; 

for  (i=l;  i<=N;  i++)  /*  Initializes  M  matrices  to  zero  */ 

{ 
for(j=l;j<=3;j++) 

{ 

M[i]£j].r=0.0;  M[i][j].c=0.0; 
Ms[i][j].r=0.0;  Ms[i][j].c=0.0; 
}  /*  End  j  loop*/ 

}  /*  End  i  loop  */ 

for(i=l;  i<=N;  i++) 

{ 

V[i].r=0.5*fabs((cos(2.0*i*dx)))-(e*F*i*dx); 

V[i].c=0.0; 

} 
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fcr  (i=l;  i<=N;  i++)  /*  This  initializes  the  M  matrices  in  sparse  format  */ 

{ 
if(i=l) 

{ 
M[i][2].r=1.0; 

M[i][2].c=-1.0*(cl*V[i].r+2*cl*c2); 

Ms[i][2]=Cconj(M[i][2]); 

M[i][3].r=0.0; 

M[i][3].c=cl*c2; 

Ms[i][3]=Cconj(M[i][3]); 
}  /*  End  if  */ 

else  if  (i==N) 
{ 

M[N][l].r=0.0; 

M[N][l].c=cl*c2; 

Ms[N][l]=Cconj(M[N][l]); 

M[N][2].r=1.0; 

M[N][2].c=-1.0*(cl*V[i].r+2*cl*c2); 

Ms[N][2]=Cconj(M[N][2]); 
}  /*  End  else  if*/ 

else 

{ 
M[i][l].r=0.0; 

M[i][3].r=0.0; 

M[i][l].c=cl*c2; 
M[i][3].c=cl*c2; 
Ms[i][l]=Cconj(M[i][l]); 
Ms[i][3]=Cconj(M[i][3]); 
M[i][2].r=1.0; 

M[i][21.c=-l*(cl*V[i].r+2*cl*c2); 
Ms[i][2]=Cconj(M[i][2]); 
}  /*  End  else  */ 

}  /*  End  i  loop*/ 


for  (i=l;  i<=N;  i++) 

{ 

psi[i].r=0.0; 

psi[i].c=0.0; 


/*  End  i  loop  */ 
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for  (1=1;  i<=N;  i++) 

{ 
psi[i].r=0.6*exp(-l*(i-5001.0)*(i-5001.0)/(2*sigma)); 

psi[i].c=0.0;  /*  Remember  -  make  divisor  float  */ 

psit[i]=Ceq(psi[i]); 

}  /*  End  for  i  */ 

for(j=l;j<=N;j++) 

{ 

conj[j]=Cconj(psi[j]); 

magipsi[j]=Cmult(conj[j],psi[j]); 

}  /*  End  j  loop  */ 


/*LU  Decomp  procedure  */ 

for  (i=l;  i<=N;  i++) 

{ 
L[i][l].r=1.0; 

L[i][l].c=0.0;  /*  Initializes  Lower  main  diagonal  */ 

}  /*  End  fori  */ 

U[l][l]=Ceq(Ms[l][2]);       /*  Initializes  first  Upper  value  */ 

L[l][2]=Cdiv(one,U[l][l]);       /*  Initializes  first  lower  value  */ 
L[l][2]=Cmult(L[l][2],Ms[2][l]); 


for(j=l;j<=N;j++) 

{ 

U[j][2]=Ceq(Ms[JP]);   /*  Solves  next  2  Upper  values  */ 

temp=Cmult(L[j][2],U[j][2]); 

U[j+I][l]=Csub(Ms[j+l][2],temp); 


L[j+l][2]=Cdiv(one,U[j+l][l]);   /*  Calculates  next  Lower  value  */ 
L[j+l][2]=Cmult(L[j+l][2],Ms[j+l][l]); 
}  /*  End  j  loop  */ 
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for  (loop=0;  loop<= 16000;  loop++) 

{ 

for  (i=l;  i<=N;  i++) 

{ 

sum.r=0.0;  sum.c=0.0;         /*  This  creates  "b"  vector  */ 

if(i=l) 

{ 

temp=Cmult(M[i]  [2]  ,psit[  1  ]); 

sum=Cadd(sum,temp); 

templ=Cmult(M[i][3],psit[21); 

sum=Cadd(sum,templ); 

}  /*  End  if  =7 

else  if  (i==N) 

{ 
temp=Cmult(M[i][l],psit[N-l]); 

sum=Cadd(sum,temp); 

templ=Cmult(M[i][2],psit[N]); 

sum=Cadd(sum,templ); 

}  /*  End  else  if  */ 

else 

{ 
for(j=l;j<=3;j++) 

{ 

temp=Cmult(M[i][j],psit[(i-l)+(j-l)]); 
sum=Cadd(sum,temp); 
}  /*Endj*/ 

}  '  /*  End  else  */ 

b[i]=Ceq(sum); 
}  /*  End  i  */ 


y[l]=Cdiv(b[l],L[l][l]); 

for  (i=2;i<=N;  i++) 

{ 

sum.r=0.0;  sum.c=0.0; 


temp=Cmult(L[i-l][2],y[i-l]);  /*  Solves  intermediate  matrix  */ 
sum=Csub(b[i],temp); 
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temp=Cdiv(one,L[i]  [  1  ]); 

y[i]=Cmult(temp,sum); 

}  /*  End  i  loop  */ 


npsi[N]=Cdiv(y[N],U[N][l]); 

for  (i=N-l;i>=l;i--) 

{ 

sum.r=0.0;  sum.c=0.0; 

temp=Cmult(U[i][2],npsi[i+l]);  /*  Solves  final  answer  */ 
sum=Csub(y[i],temp); 

temp=Cdiv(one,U[i][l]); 

npsi[i]=Cmult(temp,sum); 

}  /*  End  i  loop  */ 


for(j=l;j<=N;j++) 
{ 

nconj[j]=Cconj(npsi[j]); 
nmagpsi[j]=Cmult(nconj[j],npsi[j]); 
}  /*  End  j  loop  */ 


if(loop%100==0) 

{ 
for  (j=l;  j<=N;  j++) 

{ 

printf("%.3f  ",nmagpsi[j].r); 

} 
printf("  \n"); 


} 


for(j=l;j<=N;j++) 

/ 

psit[j]=Ceq(npsi[j]); 

}  /*  End  j  loop  */ 
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}  /*  End  loop  loop  */ 

return; 

}  /*  End  program  */ 

struct  complex  Cadd  (struct  complex  x,  struct  complex  y) 

{ 

struct  complex  z; 


z.c=x.c+y.c; 
z.r=x.r+y.r; 
return  z; 

} 

struct  complex  Csub  (struct  complex  x,  struct  complex  y) 

{ 

struct  complex  z; 

z.c=x.c-y.c; 
z.r=x.r-y.r; 
return  z; 

} 

struct  complex  Cmult  (struct  complex  x,  struct  complex  y) 

{ 

struct  complex  z; 


z.c=x.r*y.c  +  y.r*x.c; 
z.r=x.r*y.r  -  x.c*y.c; 


return  z; 

\ 

j 

struct  complex  Cdiv  (struct  complex  x,  struct  complex  y) 

{ 

struct  complex  z; 

double  d; 
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d=(1.0/((y.r*y.r)+(y.c*y.c))); 

z.r=d*((x.r*y.r)+(x.c*y.c)); 

z.c=d*((x.c*y.r)-(x.r*y.c)); 

return  z; 
} 

struct  complex  Ceq  (struct  complex  x) 


{ 

struct  complex  z; 

z.r=x.r; 

z.c=x.c; 

return  z; 

} 

struct  complex  Cconj  (struct  complex  x) 

{ 

struct  complex  z; 

z.r=x.r; 

z.c=-1.0*x.c; 

return  z; 
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